/* Copyright (C) 2010 LinBox
 * Written by <brice.boyer@imag.fr>
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
 * Boston, MA 02110-1301, USA.
 */

/** \file examples/matmul-int.C
  \brief NO DOC
  \ingroup examples
  */
//#include "linbox-config.h"

#include <iostream>
#include <omp.h>

//#include <givaro/givtimer.h>

//#include "linbox/field/modular.h"
// #include "linbox/field/modular-double.h"
#include "linbox/field/modular-balanced-double.h"
#include "linbox/field/gf2.h"
#include "linbox/blackbox/dense.h"
#include "linbox/solutions/rank.h"
#include "linbox/util/matrix-stream.h"
#include "linbox/algorithms/blas-domain.h"
#include <linbox/field/givaro-zpz.h>
#include "linbox/algorithms/cra-domain.h"
#include "linbox/algorithms/matrix-hom.h"
// #include "linbox/integer.h"

#include "linbox/algorithms/cra-full-multip-fixed.h"


#include "util/timer.h"
#include "sys/time.h"

using LinBox::Integer ;

#define _LB_MM_TIMING
//#define _LB_MM_READING
#define _LB_LOG2 0.6931471807
// #define DEBUG
// #define _LB_DEBUG
#define _LB_REPET 5

typedef LinBox::ModularBalanced<double> Modd;
// typedef LinBox::Modular<double> Modd;

template<typename Mat>
Integer& MatMax(Integer& m, const Mat& A)
{
	m = 0;
	for(typename Mat::ConstRawIterator it = A.rawBegin(); it != A.rawEnd(); ++it)
	{
		if (*it > m)
			m = *it;
		else if (*it < -m)
			m = -(*it);
	}
	return m;
}



size_t IgnobleM;
size_t IgnobleN;

struct FixedBlasMatrixDouble : public LinBox::BlasMatrix<double>
{
	FixedBlasMatrixDouble() :
		LinBox::BlasMatrix<double>(IgnobleM, IgnobleN)
	{}

	template<class T1, class T2>
	FixedBlasMatrixDouble(const T1& t1, const T2& t2) :
		LinBox::BlasMatrix<double>(t1,t2)
	{}
};

#ifdef _LB_MM_TIMING
#ifdef _OPENMP
typedef LinBox::OMPTimer Mytime;
#else
typedef Givaro::Timer Mytime;
#endif
#endif
using namespace LinBox;
using namespace std;

typedef BlasMatrix<Integer> IntegerBlasMatrix ;

/// rank or rank mod p
int test (int argc, char **argv)
{
	//     commentator.setMaxDetailLevel (-1);
	//     commentator.setMaxDepth (-1);
	//     commentator.setReportStream (std::cout);


	PID_integer ZZ;


#ifdef _LB_MM_READING
	std::ifstream inputA (argv[1]);
	if (!inputA) { cout << "Error opening matrix file " << argv[1] << endl; return -1; }

	std::ifstream inputB(argv[2]);
	if (!inputB) { cout << "Error opening matrix file " << argv[2] << endl; return -1; }

	size_t rowsA, colsA, elemA, rowsB, colsB, elemB;
	//     Integer q(argv[3]);
	//     Field F(q);
#ifdef _LB_MM_TIMING
	Mytime initch; initch.clear(); initch.start();
#endif
	inputA >> rowsA >> colsA; elemA = rowsA*colsA;
	IntegerBlasMatrix A(rowsA,colsA);
#ifdef _LB_MM_TIMING
	initch.stop();
#endif
	cout << "A is " << A.rowdim() << " by " << A.coldim() ;

#ifdef _LB_MM_TIMING
	cout << ' ' << initch ;
#endif
	cout << endl;

#ifdef _LB_MM_TIMING
	initch.clear(); initch.start();
#endif
	inputB >> rowsB >> colsB; elemB = rowsB*colsB;
	IntegerBlasMatrix B(rowsB,colsB);
	//MatrixStream< PID_integer > msB ( ZZ, inputB );
	//IntegerBlasMatrix B (msB);
#ifdef _LB_MM_TIMING
	initch.stop();
#endif
	cout << "B is " << B.rowdim() << " by " << B.coldim() ;
#ifdef _LB_MM_TIMING
	cout << ' ' << initch ;
#endif
	cout << endl;
	// #pragma omp parallel sections
	{
		// #pragma omp section
		{
#ifdef _LB_MM_TIMING
			Mytime initchA; initchA.clear(); initchA.start();
#endif
			Integer * iter=A.getWritePointer();
			for(size_t i=0; i<elemA; ++i, ++iter)
				inputA >> *iter;
			//     MatrixStream< PID_integer > msA ( ZZ, inputA );
			//     IntegerBlasMatrix A (msA);
#ifdef _LB_MM_TIMING
			initchA.stop();
			cout << "A completed in " << initchA << endl;
#endif
			// A.write(std::cout<< "A:=" , ZZ, true) << ';' << std::endl;
		}

		// #pragma omp section
		{
#ifdef _LB_MM_TIMING
			Mytime initchB; initchB.clear(); initchB.start();
#endif
			Integer * iter=B.getWritePointer();
			for(size_t i=0; i<elemB; ++i, ++iter)
				inputB >> *iter;
#ifdef _LB_MM_TIMING
			initchB.stop();
			cout << "B completed in " << initchB << endl;
#endif
			// B.write(std::cout << "B:=", ZZ,true) << ';' << std::endl;
		}
	}

#else // on donne juste des tailles
	size_t dm = atoi(argv[1]); // <dm lignes/cols
	size_t mx = atoi(argv[2]); // <2^mx
	size_t colsA = 1+ dm ; //Integer::random((Integer)dm);
	size_t rowsA = 1+ dm ; // Integer::random((Integer)dm);
	size_t colsB = 1+ dm ; // Integer::random((Integer)dm);
	IntegerBlasMatrix A(rowsA,colsA);
	size_t elemA = rowsA * colsA ;
	IntegerBlasMatrix B(colsA,colsB);
	size_t elemB = colsA * colsB ;
	Integer * iter=A.getWritePointer();
	for(size_t i=0; i< elemA; ++i, ++iter){
		*iter = Integer::random(mx);
		//                assert(!std::isnan((int)*iter));
		//                std::cout << *iter << " " ;
	}
	std::cout << std::endl;
	iter = B.getWritePointer() ;
	for(size_t i=0; i< elemB; ++i, ++iter) {
		*iter = Integer::random(mx);
		//                assert(!std::isnan((int)*iter));
		//                std::cout << *iter << " " ;
	}
	std::cout << std::endl;
	cout << "A is " << A.rowdim() << " by " << A.coldim() << endl;
	cout << "B is " << B.rowdim() << " by " << B.coldim() << endl;
#endif

#ifdef _LB_MM_TIMING
	Integer mA ;
	MatMax(mA,A);
	Integer mB ;
	MatMax(mB,B);

	long double cmplxt = (elemA * B.coldim()*std::max(logtwo(mA),logtwo(mB)))/1e6;
#endif


	for (size_t R = 0 ; R < 2 ; ++R)
	{
		{
			size_t m = A.rowdim();
			size_t n = B.coldim();
			if (!A.coldim()) m = n = 0 ;
			IntegerBlasMatrix D(m,n);
#ifdef _LB_MM_TIMING
			Mytime crachr; crachr.clear(); crachr.start();
#endif

			CRAIntegerMul(D,A,B);
#ifdef _LB_MM_TIMING
			crachr.stop();

			std::cout << "\033[0;36m * \033[0mmatmul with chinese remaindering: (" << crachr << ")\033[0;35m " << cmplxt /crachr.usertime() << "\033[0m Mflops" <<std::endl;
#endif

			bool ok = CheckMatMul(ZZ, D,A,B);
			if (!ok) std::cout << "\033[0;31m *  \033[0mproduit faux !" << std::endl;
			else std::cout << "\033[0;32m * \033[0mproduit ok !" << std::endl;

		}
		{
			IntegerBlasMatrix D(A.rowdim(),B.coldim());
			MatrixDomain<PID_integer> BMD(ZZ);

#ifdef _LB_MM_TIMING
			Mytime chrono; chrono.clear(); chrono.start();
#endif
			for (size_t r = 0 ; r < _LB_REPET ; ++r)
				MyIntegerMul(D, A, B);
#ifdef _LB_MM_TIMING
			chrono.stop();
			std::cout << "\033[0;36m * \033[0mhome cooked matmul : (" << chrono << ")\033[0;35m " << cmplxt/chrono.usertime()*_LB_REPET << "\033[0m Mflops" <<std::endl;
#endif
			cout << "D is " << D.rowdim() << " by " << D.coldim() << std::endl;
			// D.write(std::cout << "D:=", ZZ, true) << ';' << std::endl;

			Integer mD; MatMax(mD, D);
			std::cout << "D max: " << logtwo(mD) <<  " (" << _LB_LOG2*logtwo(mD) << ')' << std::endl;

			bool ok = CheckMatMul(ZZ,D,A,B);
			if (!ok) std::cout << "\033[0;31m *  \033[0mproduit faux !" << std::endl;
			else std::cout << "\033[0;32m * \033[0mproduit ok !" << std::endl;

			// {
			// for(IntegerBlasMatrix::RawIterator it = D.rawBegin(); it != D.rawEnd(); ++it)
			// std::cout << *it << " ";
			// std::cout << std::endl;
			// }

		}
		{
			IntegerBlasMatrix D(A.rowdim(),B.coldim());
			MatrixDomain<PID_integer> BMD(ZZ);

#ifdef _LB_MM_TIMING
			Mytime chrono; chrono.clear(); chrono.start();
#endif
			for (size_t r = 0 ; r < _LB_REPET ; ++r)
				BMD.mul(D, A, B);
#ifdef _LB_MM_TIMING
			chrono.stop();
			std::cout << "\033[0;36m * \033[0mmatmul over integer : (" << chrono << ")\033[0;35m " << cmplxt/chrono.usertime()*_LB_REPET << "\033[0m Mflops" <<std::endl;
#endif
			cout << "D is " << D.rowdim() << " by " << D.coldim() << std::endl;
			// D.write(std::cout << "D:=", ZZ, true) << ';' << std::endl;

			Integer mD; MatMax(mD, D);
			std::cout << "D max: " << logtwo(mD) <<  " (" << _LB_LOG2*logtwo(mD) << ')' << std::endl;

			bool ok = CheckMatMul(ZZ,D,A,B);
			if (!ok) std::cout << "\033[0;31m *  \033[0mproduit faux !" << std::endl;
			else std::cout << "\033[0;32m * \033[0mproduit ok !" << std::endl;


		}
	}
	return 0;
}

int main (int argc, char **argv)
{
	if (argc != 3) {std::cout << "bad arguments" << std::endl; exit(-1) ;}
	Integer::seeding();
	for (size_t r = 0 ; r < 3 ; ++r)
		test(argc,argv);
}

#undef _LB_MM_TIMING
#undef _LB_MM_READING
#undef _LB_LOG2
#undef _LB_DEBUG
#undef _LB_REPET
/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s:syntax=cpp.doxygen:foldmethod=syntax
